Relational data and table joins

• Much of the world’s data lives in relational database management systems (RDBMSes), or simply “relational databases”
• Sometimes you might work with one of these directly; other times you might receive CSV files which originated in a RDBMS
• “Relational” == Table A has a relationship to Table B, via column X which appears in both tables (“foreign key relationship”)
• Why relational?
  • One big table gets to be unmanageable, but we also don’t want to duplicate info in the multiple places that we need it
  • Often-repeated values (e.g., a drug name) can be more efficiently stored/retrieved if they are replaced by a number, which points to the human-readable name.

Joining synthesizes a new a table from two or more parent tables.

Relationship types

• One-to-one relationships: one row in Table A corresponds to one row in Table B
• One-to-many relationships: “lookup” function; controlled vocabulary; Table B describes categories that are referenced in Table A
• Many-to-many relationships: there are many subjects (Table A), and many genetic markers (Table B); a patient may have an unspecified number of markers and vice versa, a relationship described in Table C

It can be helpful to have these relationships in mind when thinking about different ways we might want to join tables together.

Types of joins

Join type asks, “if I have unpaired rows when I join Tables A and B, what do I do with those rows?”

• Full join: keep all rows from both tables, filling NA as needed
• Inner join: keep only rows that appear in both tables
• Left join: inner join + keep all of the rows from the left table, filling NA as needed
• Right join: inner join + keep all of the rows from the right table, filling NA as needed
• Filtering joins: filter rows of Table A based on presence/absence of the row’s ID in Table B

Each of the above joins might produce a different number of rows.

semi_join(), anti_join(): filter x vs. y

Filtering joins decide which rows to keep from x, based on the presence (or absence) of the value in y

Example: you applied inclusion criteria to get a list of patient ID’s to include in analysis, and now you would like to filter just those patients’ records into a data frame.

inclusion_ids <- tibble(subject_id = c(9134, 663, 5408))

covid_testing %>% 
  semi_join(y = inclusion_ids, by = join_by(subject_id))

# A tibble: 3 × 17
  subject_id fake_…¹ fake_…² gender pan_day test_id clini…³ result demo_…⁴   age
       <dbl> <chr>   <chr>   <chr>    <dbl> <chr>   <chr>   <chr>  <chr>   <dbl>
1       9134 grunt   rivers  male         7 covid   clinic… negat… patient   0.8
2        663 ithoke  targar… male         9 covid   clinic… negat… patient   0.8
3       5408 alia    ryswell female      10 covid   clinic… negat… patient   0.9
# … with 7 more variables: drive_thru_ind <dbl>, ct_result <dbl>,
#   orderset <dbl>, payor_group <chr>, patient_class <chr>, col_rec_tat <dbl>,
#   rec_ver_tat <dbl>, and abbreviated variable names ¹​fake_first_name,
#   ²​fake_last_name, ³​clinic_name, ⁴​demo_group

Jobs involving factors

Typical things you’ll need to do with factors:

• Convert character vector into factor: as_factor()
• Change the label associated with a level: fct_recode()
• Collapse two or more levels into one level (e.g., “Other”): fct_collapse()
• Change factor’s order of levels: several functions

my_fct <- as_factor(c("A", "B", "C", "B", "A", "B"))

my_fct

[1] A B C B A B
Levels: A B C

my_fct %>% fct_recode(beta = "B")

[1] A    beta C    beta A    beta
Levels: A beta C

my_fct %>% fct_collapse(Other = c("A", "B"))

[1] Other Other C     Other Other Other
Levels: Other C

my_fct %>% fct_infreq()

[1] A B C B A B
Levels: B A C

Pivoting: what and why

• To “pivot” tabular data means to restructure part of it in terms of column and row definition
• Example: a variable described by one column per year of interest
  • But ggplot expects all values to be in the same column
  • Solution: all values into a new value column, and each year into a new year column
  • \(n \mathrm{\ columns} \rightarrow 2 \mathrm{\ columns}\)
• The goal is often “tidy data”: one observation per row, one variable per column

pivot_longer()

• For when you need reduce the number of columns, making the table longer with more rows
• Or: when multiple observations are compressed into a single row
• Note required arguments (given in example below)
• Example: annual teen birth counts by state, one column for each year (2013-19); source: CDC
  • cols = the year columns (`2013`:`2019`)
  • names_to =new column where you want the years recorded (year)
  • values_to =new column where you want the values recorded (births)

Checking for normality: histogram + theoretical

The ggplot2 stat_function() layer allows you to plot various statistical functions (such as a theoretical distribution). Note that for the two layers to be compatible, we transform the histogram’s y to density. dnorm() is the base R function for the Normal distribution.

mu <- mean(polyps$number12m, na.rm=TRUE)
sigma <- sd(polyps$number12m, na.rm=TRUE)

polyps %>% 
  drop_na(number12m) %>% 
  ggplot() +
  geom_histogram(aes(x=number12m, y=after_stat(density)), binwidth = 3) +
  stat_function(fun=dnorm, geom="line", color="blue", linewidth=1.5,
                args=list(mean=mu, sd=sigma))

Inference and modeling - the base R way

• A variety of functions are available depending on the test / probability distribution / regression you need
• These functions usually expect vectors as input; usually produce some printable object (summary() also usable)
• Student’s t-test: t.test()
• Chi-squared test: chisq.test() (note: expects a matrix for 2-way)
• Linear model: lm() with predict()

These functions can be compatible with tidyverse work (but not %>% )

Are mean 12-month counts significantly different between sulindac (treatment) and placebo groups?

sulindac12m <- polyps %>% 
  filter(treatment == "sulindac") %>% 
  pull(number12m)
sulindac12m

 [1] NA  2 17  1 25  3 33 NA  3  1  4

placebo12m <- polyps %>% 
  filter(treatment == "placebo") %>% 
  pull(number12m)
placebo12m

 [1] 63 28 61  7 15 44 28 10 40 46 50

t.test(sulindac12m, placebo12m)

    Welch Two Sample t-test

data:  sulindac12m and placebo12m
t = -3.6114, df = 16.901, p-value = 0.002172
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -40.79597 -10.69898
sample estimates:
mean of x mean of y 
 9.888889 35.636364 

Is there a relationship between baseline and 12-month counts?

polyps_model <- lm(number12m ~ baseline, data=polyps)
polyps_model

Call:
lm(formula = number12m ~ baseline, data = polyps)

Coefficients:
(Intercept)     baseline  
    19.1647       0.1113  

plot(polyps$baseline, polyps$number12m)
abline(polyps_model, col="blue")

What 12-month count does this model predict for a baseline of 50? of 100?

new_baseline <- tibble(
  baseline = c(50, 100)
)
predict(polyps_model, newdata=new_baseline)

       1        2 
24.72882 30.29290 

Inference and modeling - the tidymodels way

• Specifically, the {infer} and {parsnip} packages
• New verbs for the tidyverse pipeline: (or you can use wrappers like t_test())
  • specify() variable/relationship of interest
  • hypothesize() the null hypothesis
  • generate() (or simulate if theoretical) the null distribution
  • calculate() a test statistic like Chisq, F, t, or z
• Good vignettes for getting started: Getting to Know infer, Introduction to parsnip

Are mean 12-month counts significantly different between sulindac (treatment) and placebo groups?

# infer::t_test()

polyps %>% 
  t_test(response = number12m,
         explanatory = treatment)

# A tibble: 1 × 7
  statistic  t_df p_value alternative estimate lower_ci upper_ci
      <dbl> <dbl>   <dbl> <chr>          <dbl>    <dbl>    <dbl>
1      3.61  16.9 0.00217 two.sided       25.7     10.7     40.8

Is there a relationship between baseline and 12-month counts?

# parsnip::linear_reg()

polyps_model <- linear_reg() %>% 
  set_engine("lm") %>% 
  fit(number12m ~ baseline, data = polyps)
polyps_model

parsnip model object


Call:
stats::lm(formula = number12m ~ baseline, data = data)

Coefficients:
(Intercept)     baseline  
    19.1647       0.1113  

What 12-month count does this model predict for a baseline of 50? of 100?

new_baseline <- tibble(
  baseline = c(50, 100)
)
predict(polyps_model, new_baseline)

# A tibble: 2 × 1
  .pred
  <dbl>
1  24.7
2  30.3

Plotting a model (without geom_smooth)

Run predict() on a data frame of your explanatory values. Then bind the predicted column back with the original data frame and plot the predictions with their own geom_line().

baseline_df <- polyps %>%
  select(baseline)

predicted_df <- predict(polyps_model, new_data = baseline_df) %>% 
  rename(predicted_12m = ".pred") %>% 
  mutate(predicted_12m = round(predicted_12m, 1))
predicted_df

# A tibble: 22 × 1
   predicted_12m
           <dbl>
 1          19.9
 2          27.7
 3          19.9
 4          19.7
 5          21.7
 6          23.1
 7          20.4
 8          20.5
 9          19.9
10          54.6
# … with 12 more rows